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Abstract 

An  implementation,  using  Gaussian  LU  decomposition 
with  row  interchanges, of  Stiefel's  exchange  algo¬ 
rithm  for  determining  a  Chebyshev  solution  to  an 
overdetermined  system  of  linear  equations  is  pre¬ 
sented.  The  implementation  is  computationally 
more  stable  than  those  usually  given  in  the  lit¬ 
erature.  A  generalization  of  Stiefel's  algorithm 
is  developed  which  permits  the  occasional  exchange 
of  two  equations  simultaneously.  Finally,  some 
experimental  comparisons  are  offered. 
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1.  Introduction 


The  problem  of  finding  a  vector  x  ■  (x^,...,  xn)  which  solves  an 
overdetermined  system  of  equations 

n 

x)  3  E  a 
d-i 

in  the  sense  that 


ijXj 


-  d. 


(i=l,...,  m;  where  m  >  n) 


max 

i.<i<m 


max 

l<i<m 


|rt(x)| 


for  any  x  e  En  is  treated  by  Stiefel  in  [1].  Such  an  x  is  called  a 
Chebyahev  or  mlnlmax  solution  to  the  system. 

Given  an  overdeterrained  system  of  linear  equations  Ax  =  d  whose 
matrix  of  coefficients  satisfies  the  Haar  condition  (each  n  X  n  sub¬ 
matrix  is  nonsingular),  Stiefel  presents  in  [l]  an  algorithm  called  the 
exchange  method  for  finding  a  Chebyshev  solution.  In  a  later  paper, 

[2],  the  exchange  method  is  shown  to  be  equivalent  to  the  simplex  method 
applied  to  a  suitable  linear  programming  problem. 

In  this  regard,  Stiefel  suggests  the  use  of  techniques  drawn  from 
the  simplex  method  for  the  implementation  of  his  algorithm.  These 
techniques  are  characterized  by  their  use  of  Jordan  elimination,  for  the 
most  part  without  row  or  column  interchanges  to  pick  the  most  advanta¬ 
geous  pivots,  for  solving  linear  equation  systems  which  arise  during  the 
computation.  These  methods  are  fast  but  computationally  unstable.  In 
this  paper  we  propose  a  computational  scheme  based  upon  the  more  stable 

1 


method  of  Gaussian  LU-decomposition  using  row  interchanges.  Attention 
is  paid  to  the  peculiarities  of  the  exchange  method  to  make  computation 
as  fast  as  possible. 

Afterwards  a  generalization  of  Stiefel's  algorithm  is  presented 
which  permits  the  occasional  exchange  of  two  equations  at  once. 

Finally  some  experimental  comparisons  of  selection  rules  for  use 
with  the  exchange  method  are  tabulated. 


2.  Background  Theory 


There  is  a  full  treatment  of  the  theory  and  the  exchange  method  in 
Chapter  2  of  [9].  (The  exchange  method  is  called  the  ascent  algorithm 
in  this  work. )  We  therefore  confine  ourselves  in  this  section  and  the 
next  to  a  statement  of  pertinent  results,  omitting  proofs. 

According  to  corollary  7.4.7.,  page  410,  of  [4],  any  overdetermined 
system  of  linear  equations  has  a  Chebyshev  solution.  The  following 
lemma  and  theorem  serve  to  characterize  these  solutions. 


Lemma: 


Let  B  =  [b  ]  be  a  p  X  q  matrix  with  rows  B-*...,  B 
1  j  ip 

There  is  a  vector  y  *  (y^,...,  y^)  such  that 


<  0 


for  all  i-1, . . . ,  p 


P 

if  and  only  if  0  ^  T  for  all  nontrivial  choices  of 

f»l 

of  nonnegative  scalars  G^, . . . ,  0^  . 

This  lemma  is  a  special  case  of  corollary  6,  page  115,  of  [5]. 
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Let  Ax  ■  d  be  an  overdetermined  system  of  m  linear  equations  in 
n  unknowns.  For  any  vector  x  ■  (x^,...,  xr),  denote  the  residuals 

(1-1...-,  »)  by  r^x)  . 

4*  h 

Let  A i  be  the  i  row  of  the  matrix  A  . 

Given  any  fixed  vector,  v  »  (v^,...,  v^),  we  may  assume  with  no  loss 
of  generality  that  the  equations  have  been  ordered  and  numbered  so  that 

l<Km  lri(v)l  ‘  lrl<T)l  "  •"  *  lrk(v)l  >  lrkn<v)l  i  •••  -  |rn(v)  I  ’ 

where  1  <  k  <  m  . 

Theorem:  There  is  a  vector  z  for  which 

X<l<n  K«l  < 
k 

if  and  only  if  0  /  £  w^  sgn(r^(v))A^  for 

1*1 

all  nontrivial  choices  of  nonnegative  scalars 

wL, •  •  • ,  wR  • 

For  the  purposes  of  the  exchange  method  we  restrict  our  attention 
henceforth  to  overdetermined  systems  of  m  linear  equations  in  n  un¬ 
knowns,  Ax  ■  d,  for  which  rank(A)  =  n  . 

To  begin,  suppose  that  m  =  n+l  .  There  is  no  loss  of  generality  in 
assuming  that  the  equations  have  been  ordered  so  that  the  first  n  rows, 

A. ,  •  •  • ,  A  ,  of  A  are  linearly  independent.  Thus,  scalars  \ 

in  i  n+l 

can  be  found  with  \  . ,  4  0  such  that 

n+l  ' 


n+1  n+1 

And  so  0  »  sgn(e)  0  =  £  l^i  I si  s6n(«)Ai  ■  ^  |X±  |sgn(ri(x))Ai  . 

Hence,  by  the  preceding  theorem,  x  =  xr)  is  a  Chebyshev  solu¬ 

tion  for  the  given  system.  (For  an  alternate  discussion  of  (n+l)  X  n 
systems  see  [6].) 

Returning  to  the  general  case  (m  >  n+l),  suppose  for  some  set  of 
n+1  rows  of  A  the  first  n  of  which  are  linearly  independent  (with 
complete  generality,  the  first  n+1  rows  of  A)  we  construct  the 


k 


Following  the  convention  put  forth  in  [l],  any  subsystem 
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of  the  given  system  with  rank 


uur 


will  called  a  ref¬ 


erence  .subsystem,  and  the  rows  A.  >...,  A.  will  be  called  a  refer- 

1  n+1 

ence  set.  If  x  ■  (x^,...,  xn)  is  a  Chebyshev  solution  to  a  reference 
•.*  'bsysteiUj.  the  value. 

M  *  |rt  (x) -  |r  (x) |  -  inf  K^+1lri  Ml 
X1  n+1  ycEn  -J-ni  xj 


will  oe  called  the  reference  deviation  for  the  reference  subsystem.  It 
i,  '>'l.a.quely  determined  by  the  reference  subsystem. 


3  The  Exchange  Method 

StJpfel's  algorithm  consists  of  starting  with  a  reference  subsystem 
ard  modifying  it  one  equation  at  a  time  so  as  to  increase  the  reference 
deviation  by  each  change.  Each  modification  proceeds  as  follows: 

is  a  reference  set.  Let 
x  =  (x1,...,  xn)  be  a  Chebyshev  solution  to  the  corresponding  reference 
sui^vst^m  computed  as  above.  So  we  have  t,  \^>...>  Xn+^  which  satisfy 

n+1 

a)  £  XjA.  =  0 


We  may  assume  that  . . . , 


c)  r^(x)  ■  si«  for  i  «  1,...,  n+1, 
where  s^  ■  sgn(\^)  • 

If  x  is  not  a  Chebyshev  solution  to  the  full  given  system,  then  by 
the  discussion  in  the  previous  section,  there  is  an  ««{n+2,...,  m} 
for  which  |r^(x)|  >  |c|  •  Let  pn+^  be  scalars  for  which 

n+1 

'  A«  '  &  PiAl  ' 

In  order  to  proceed,  we  impose 

Condition  1:  ^  0  for  all  i  =  1,.*.,  n+1  . 

If  this  holds,  let  9c(l,...,  n+l}  be  such  that 


gaSPB  max  V^i 
Xp  “  l<i<n+l  * 

where  o„  -  sgn(r^(x)),  and  s  =  sgn(*)  . 

are  a  reference  set. 

to  the  reference  subsystem 


Now  impose 

Condition  2;  Ap_1,  A  kn+±,  A& 

We  form  a  Chebyshev  solution  x'  *=  (x^,. . * ,  x^) 


«sTi. 

i _ 

- 1 

•  N 
(— * 

1 

- 1 

_ 1 

V 

• 

s 

Vi 

Vi 

• 

• 

^(5+1 

• 

V, 

z 

L  nJ 

Vi 

1 

— 1 

•o* 

_ 1 

in  the  usual  fashion, 
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producing  c ' ,  X^> •••> 


s-i*  s+i . xa+i<  x; such  that 


n+l 

a')  £  x;  Ai  +  x; 

Me 


'A  *  0 
Of  a 


n+l 

Z  *  Kdc 


lx;' +  l^j1 

Me 

c * )  ri(x')  *  s|«'  for  1*1,. 
where  =  sgn(x|)  • 


8+l)«  •  •  )  11+1)  Qt 


We  further  have 


xi  “  Vixi  [  "xi2  ’  3  n+1;  Me)  • 


(Note  that,  by  the  choice  of  0,  the  product  of  the  term  in  brackets 
with  s  *  sgn(c)  is  nonnegative.) 


Furthermore, 


if  K  = 


n+l  | X  | 

lx’ I  +  £  IxJ  and  c  -  -Lil-  , 

®  k-i  *  K 

Me 


it  can  readily  be  shown  that 


I « '  I  =  c|ra(x)|  +  (l-c)|e|  . 
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It  is  important  to  note  that>  if  condition  1  is  satisfied  by  the 

second  reference  set  (i.e.,  ^  0  for  i=l,...,  0-1,  0+1,...,  n+1,  a)> 

then  c  >  0  .  Therefore  | « '  |  >  |c|,  since  |r  (x)|  >  I  cl  •  The 

a 

strictness  of  the  inequality  |c'|  >  |«|  implies,  by  a  simple  contra¬ 
diction  argument,  that  if  an  initial  reference  set  is  chosen  and  subse¬ 
quently  modified  as  above  by  exchanging  successive  non-reference  set 
rows  of  the  matrix  A  for  rows  in  the  reference  set,  and  if  conditions 
1  and  2  hold  at  each  exchange,  the  process  must  converge  upon  a 
Chebyshev  solution  for  the  full  system. 


J+.  Jordan  Elimination 


An  excellent  example  of  an  implementation  of  the  exchange  method 
which  uses  Jordan  elimination  is  given  on  page  50  of  [9]. 

Briefly,  given  indices  fi^...,  in+1l  C  {!>•••>  m},  numbers 


•  > 


n+1 


are  found  so  that 


and 


n+1 

I 


k=l 


n+1 


XA  “ 0  • 


Setting  sk  =  sgn(\k)  for  k^l,...,  n+1,  the  matrix 


at 

\ 

•  «  t 

AT 

i 

n+1 

sl. 

•  •  • 

sn+l 
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A.  V* 

-V  -J' 


h 


-v 


is  formed  using  a  sequence  of  n+1  pivot  operations  (Jordan  elimination 
steps). 

Each  exchange  step,  then,  involves  forming 


(x^,...,  x^,  c)  =  [d^  ,••*,  d^  ]C, 

X1  n+1 


computing 


n+1 

rj  ■  Z  -  dj  for  aU  Jl11!— -  Vl’ 


selecting  a  so  that  |r  |  =  max,  and  forming 

ot 


[pl>--->  Pn+13  -  aa,n'  sS^(ra)]C 


The  last  column  of  C  has  the  form 


Vo 

X2/G 

Xn+l/G 


where 


n+1 

0  -  Z  W  * 

k=l  K 


Hence,  B  is  selected  as  an  index  for  which 


10 


max  . 


Sgn(,a)  sgn(c)  PB/0B)n+1  * 


An  appropriate  pivot  operation  on  C  ends  the  exchange  step. 
The  can  be  found  in 


4p-  +  2n 


+  1 


operations  (counting  only  multiplications  and  divisions),  and  the  initial 

3  2 

computation  of  C  requires  an  additional  n  +  Jn  +  n  operations. 

In  each  exchange  step  the  quantities 


xn+1,  e,  Pl,...,  pn+1 

2 

require  2n  +  4n  +  2  operations  to  compute,  and  the  updating  of  C 

2 

demands  an  additional  n  +  2n  +  1  operations.  Hence,  k  exchanges 
may  be  carried  out  with 


+  (3k  +  5)n2  +  (6k  +  ~)n  +  3k  +  1 


operations. 

While  row  and  column  interchanges  can  be  permitted  during  the 
initial  sequence  of  Jordan  elimination  steps  which  forms  C,  so  that 
pivot  elements  of  largest  possible  magnitude  can  be  selected,  no  pivot 
choice  is  possible  during  the  subsequent  updatings  of  C  .  For  simple 
examples  of  the  danger  implicit  in  this  fact  see  [10,11].  The  danger 
is  studied  at  greater  d-*pth  in  [3»7»83* 


11 


»  i?  vi'isi' 


I  , 

>  \ 


V 


Ji'.v 


5.  LU  Decomposition 

Starting  from  any  reference  subsystem  of  the  given  overdetermined 
system,  the  exchange  method  produces  a  new  reference  subsystem  at  the 
cost  of  solving  three  nonsingular  sets  of  n+1  linear  equations: 


<  > 

1  .  j  • 

\ 


r 


PX  =  r. 


T 

P  X  =  Tr 


Pp  =  r. 


The  vector  r^  is  given,  but  r^  depends  upon  \  and  r^  depends 
upon  x  •  If  three  such  systems  of  equations  were  given  in  isolation, 
the  general  method  of  solution  would  consist  of  making  an  accurate 
LU  decomposition  of  P  using  Gaussian  elimination  and  backsolving  six 
triangular  systems  of  linear  equations.  This  can  be  done  with 


, 

..  '< 


y-  +  4n 2  +  0(n) 


operations.  With  Stiefel's  algorithm,  however,  this  price  need  not  be 
paid  at  every  exchange.  The  matrix  P',  derived  from  P  by  one 
exchange,  differs  from  P  only  in  its  £j  column.  If  column  inter¬ 
changes  are  not  permitted  in  computing  LU  decompositions,  then  the 
decomposition,  L'U' ,  of  P'  is  identical  in  certain  portions  to  the 
decomposition,  LU,  of  P,  affording  a  saving  of  work.  Furthermore, 
pivotal  selection  using  row  interchanges  can  be  allowed.  While  an 
example  of  a  matrix  is  given  in  [7]  for  which  this  strategy  is  poor, 
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it  is  the  strategy  commonly  used  and  is  almost  always  stable  in  practice 
(e.g.,  see  comments  to  this  effect  in  [J]  and  [6]).  In  any  event 
it  is  superior  to  the  strategy  of  making  no  pivot  selection. 

The  work  done  in  carrying  out  k  exchange  steps,  involving  columns 
0^,...,  0^  of  can  be  cut  to 


(k  +  1)(t-  +  ^n2  +  7n)  - 


n  +  1 


/iii* 


7k  +  10 


operations. 

For  example,  if  0^  *  . . .  »  ,  this  becomes 

(k  +  2)^  +  +  4)n2  +  0(n), 

roughly  half  the  work  that  would  be  required  if  no  advantage  were  taken 
of  the  similarities  between  P  and  P'  . 


6.  Detailed  Outline  of  an  LU  Implementation 

1.  Select  n+1  indices  (i  in+^}  c  fl,...,  m)  so  that  the 

matrix 


Ai, 

di, 

1 

1 

• 

• 

• 

• 

• 

• 

Ai 

di 

n+1 

n+1 

15 


V 


is  nonsingular.  If  this  cannot  be  done,  terminate  with  an 
appropriate  indication.  The  user  may  then  check  whether  the 
system  Ax  =  d  can  be  sati  sfied  exactly. 


2.  Perform  the  Gaussian  reduction  of  P  into  the  product  of  a 

unit  lower  triangular  matrix  L  and  an  upper  triangular  matrix 

U  .  All  information  about  L  and  U  can  be  stored  in  the 

T 

space  initially  occupied  by  P  plus  one  vector  (for  inter¬ 
change  information).  In  each  column  the  element  of  largest 

magnitude  on  or  below  the  diagonal  is  to  be  used  as  the  pivot. 

T 

If  the  LU  decomposition  of  a  matrix  differing  from  P  only 

in  the  0  column  is  available,  one  can  save  computation  by 

using  the  first  0-1  columns  and  (as  pointed  out  by  W.  Kahan 

of  Toronto)  the  upper-right-hand  (0-l)  X  (n-0)  submatrix  of 

this  decomposition  as  the  corresponding  segments  of  the  decom- 
T  T 

position  of  P  .  If  rank(P  )  <  n+1,  terminate. 

3*  Solve 


PT 

• 

• 

=  LU 

M 

• 

85 

"  0  ' 

• 

• 

• 

Xn+1 

Xn+1 

0 

-1  • 

m  — 

This  requires  the  forward- solution  of 


followed  by  the  back- solution  of  U\  =  v  .  (Permutations 

due  to  the  row  interchanges  of  step  (2)  are  ignored  in  the 

remainder  of  the  outline).  If  v^, ...,  v^  are  available  from 

a  forward- solution  involving  an  L  whose  first  3-1  columns 

are  identical  with  those  of  the  matrix  L  being  used  here, 

only  v-,...,  v  need  be  computed.  If  any  is  zero, 

p  n  1 

terminate . 

4.  Set 


n+1 

« - 1/ 1 

i=l 


IxJ 


If  c  is  less  than  any  value  of  e  previously  computed  for 
the  current  data,  go  to  step  (9)* 

5.  Solve 


Px 


c 


sgn^) 

•  | 

sen(\,+l>  • 


x  will  turn  out  to  be  -1  . 
n+1 


6.  Compute 


r 


J 


w 


I,  ajA  -  dj 


15 


./ 


\ 


for  each 


in+l^ 


Let  0/  be  an  index  for  which  |r^(x)|  is  iraximal.  If 

|r  (x)|  <  then  (x.  x  )  is  a  candidate  as  the 

ot  —  i.  n 

Chebyshev  solution  of  Ax  =  d;  go  to  step  (10). 

7.  Solve  Pjj[  =  (oth  column  of  A1)  . 

8.  Find  0e{lj***i  n+1}  so  that 

^  E8n[ro(x)1 

is  maximal.  Replace  the  set  of  indices  fi^>...>  }  by 

tii,,,‘'  S-i*  *»  Vi'***'  Vl5  ' 


•  th  T  T 

Replace  the  0  column  of  P  by  A  .  Go  to  step  (2). 

Of 

9-  Restore  the  preceeding  set  of  indices  {i^,...,  ^n+]_}  an<^ 
recover  the  preceeding  LU  decomposition. 

10.  Iteratively  refine  the  solution  to  the  system 


PX  =  € 


sgn(\L) 

sgn(xn+i) 


16 


s _ _ 


.  > 


according  to  the  scheme  given  on  page  121  of  £33  -  (The  con¬ 
vergence  of  this  refinement  process  is  established  in  [12]). 


Check  the  residuals  *\j(x)  for 


|ra(x)l  =7  |rJ(x)l  *  e' 


then  give  [x^>...,  xn+1J  as  the  Chebyshev  solution.  If  this 
residual  check  is  not  successful,  but  the  refinement  process 
has  been  carried  out  before  and  the  last  refined  value  of  « 


is  greater  than  the  current  refined  value  of  e,  return  the 


last  refined  values  of  x, ,...,  x  as  a  doubtful  solution. 

1  n 

Otherwise  return  to  step  (7)- 


7*  Remarks  on  the  Outline 


We  have  ignored  scaling  strategies  in  programming  our  implementation. 

Step  (10)  serves  to  improve  the  final  values  of  e>  x^,...,  . 

It  is  usually  performed  only  once.  It  is  not  uncommon  to  produce  values 
for  tj  x^,...,  xn  which  are  correct  substantially  to  full  machine 
precision;  i.e.,  compare  runs  A  and  D  in  the  appendix.  The  decisions 
made  in  step  (10),  after  the  refinement,  have  been  included  as  an  attempt 
to  supply  the  Chebyshev  solution  for  the  reference  subsystem  having  the 


i  "•S’ . 


-.1  • 


largest  reference  deviation  in  those  infrequent  cases  where  the  test 


max 


l*j(x)|  <  c 


consistently  fails  to  be  satisfied. 

T 

Note  that  the  LU  decomposition  of  P  is  used  to  solve  the  system 

of  equations  Px  *=  <  sgn(x)  (step  5).  In  [3]  it  is  shown  that  the 

computed  solution  to  Ax  =  b  via  LU  decomposition  is  the  exact  solution 

to  (A  +  K)x  =  b,  where  a  bound  on  ||k||  can  be  placed.  It  is  easily 

00 

T 

shown  that  the  computed  solution  to  A  y  ■  d  via  the  LU  decomposition 
of  A  is  the  exact  solution  to  (A  +  H)y  =  d,  where  the  same  bound 
pertains  to  ||h|| 

00 

8.  Algol  60  Description 

procedure  Chebyshev  (A,d,h,m,n,refset,epz,insufficientrank,zerolambda); 
value  m,n;  integer  m,n;  real  array  A*d,h; 

Integer  array  refset;  real  epzj  label  insufficientrank,  zerolambda; 
begin 

real  procedure  ipr  (ii,ii,uu,aa,bb,cc); 
value  U,uu,ccj  real  aa,bb,cc;  Integer  ii;if,uu; 
begin  comment  single-precision  inner-product  routine; 
real  sum; 
sum  :=  cc; 

for  ii  :=  li  step  1  until  uu  do  sum  :=  Bum  +  aaxbb; 
ipr  :=  sum; 


I 

[ 


D 


t 

[ 


L 

L 

0 


end  ipr; 

real  procedure  ip2  (ii,n,uu,aa,bb,cc); 

comment  ip2  is  a  version  of  ipr  which  accumulates  the  products  aaxbb  in 
a  double-precision  sum,  whose  final  value,  rounded  to  single¬ 
precision,  is  taken  as  the  value  of  ip2.; 
procedure  trisolv  (fis,fid,fie,sis,sie,fi,si,so|,rhs,mat,piv,vip); 
value  fis,fid,fie;  integer  fis,fid,fie,sis,sie,fi,si; 
real  so!,rhs,mat,piv;  real  procedure  vip; 
begin  real  tl,t2; 

comment  trisolv  solves  a  triangular  system  of  linear  equations.  The 
off-diagonal  part  of  the  system's  matrix  is  given  by  mat,  the 
diagonal  part  by  piv,  and  the  right  hand  side  of  the  system  by 
rhs.  The  solution  is  developed  in  sol.  By  appropriately 
setting  the  first  five  parameters,  either  an  upper  or  a  lower 
triangular  system  can  be  treated.  Column-by- column  Gauss 
decomposition  of  a  matrix  can  be  compactly  expressed  using 
trisolv.  vip  is  a  vector  inner-product  routine. ; 
for  fi  :=  fis  step  fid  until  fie  do 

begin  tl  :=  -vip  (si,sis,sie,soi.,mat,-rhs ) ;  t2  :=  piv; 
si  :  =  fi;  sol  :=  _if  t2  =  1  then  tl  else  tl/t2; 


end; 

end  trisolv; 

Boolean  finished;  switch  decompbranch  :=  return, itr ; 
switch  failures  :=  insuff icientrank, zerolambda; 
integer  ml,nl,npl,i,j,k,i ,b,af  ,al,i st,/0,ll,/01,cnt; 
real  lasteps,preveps,ref ,s, t,cps,eta, cnor  n,snorm; 


real  array  P[o:n,o:n],lam,rv, sv,x,w,xr[o:n) ; 
integer  array  r[o:n],ix[o:m-l]; 

comment  The  subsystem  of  n+1  equations  currently  being  investigated 
is  listed  in  ix[o],...,  ix[n]  .  The  other  equations  are  listed 
in  the  remainder  of  ix  .  r  contains  row  indices.  Row  inter¬ 
changes  during  the  Gauss  decomposition  of  P  are  carried  out 
by  permuting  the  elements  of  r  ; 
procedure  resid  (vip);  real  procedure  vip; 
begin 

comment  resid  computes  those  components  of  the  residual  vector  Ax-d 
associated  with  the  equations  not  in  the  reference  subsystem. 
The  sign,  magnitude,  and  associated  equation  number  of  the 
largest  component  are  saved,  vip  is  a  vector  inner-product 
routine. ; 
ref  :»  -1; 

for  j  :=  npl  step  1  until  ml  do 
begin 

i  :*=  ix[j]j 

t  :»  vip  (k,0,nl,x[k],A[i,k],-d[i])j 
if  abs  (t)  >  ref  then  begi.i  ref  :=  abs  (t); 

al  :=  j;  s  :»  sign  (t); 

end; 

end; 

end  resid; 

ml  :=  m-1;  nl  :=  n-1;  npl  :«  n+1; 
lasteps  :«  0;  preveps  :=  -1; 
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for  i  :=  0  step  1  until  n  do  r[i]  :=  ix[i]  :=  i; 
for  i  :=  npl  step  1  until  ml  do  ix[i]  :=  i; 
comment  The  initial  reference  subsystem  is  chosen  by  making  a  copy  of 
the  transpose  of  A  bordered  with  d  and  carrying  out  a 
Gaussian  reduction  upon  it  with  row  and  column  interchanges 
used  to  select  the  largest  possible  pivot  at  each  stage. ; 

begin 

real  array  TAB[o:n,o:ml] ; 

for  j  :=  0  step  1  until  ml  do 
begin 

TAB[n,j]  :=  d[ j ] ; 

for  i  :=  0  step  1  until  nl  do  TAB[i,j]  :=  A[j,i]; 
end; 

for  i  :=  0  step  1  until  n  do 
begin 

t  :=  0; 

for  j  :=  i  step  1  until  n  do 
begin 

k  :=  r[j] ; 

for  I  :=  i  step  1  until  ml  do 
begin 

ref  :=  TAB[k,ix[l]]; 
if  abs  (ref)  >  t  then 

begin  s  :=  ref;  t  :=  abs  (ref);  a /  :=  j;  b  :=  I;  end; 
end; 
end; 
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if  t  =  0  then  begin  j  :=  1;  £0  to  singular;  end; 
k  :=  r[a/];  r[ai]  :=  r[i];  1st  :=  r [i ]  :=  k; 
k  :=  ix[b];  ix[b]  :=  ix[i];  al  :=  ix[i]  :=  k; 
for  j  :=  i+1  step  1  until  ml  do 
begin 

t  :=  ix[ j ] ; 

ref  :=  TAB[ist,i]/s; 

for  k  :=  i+1  step  1  until  n  do 

begin 

a I  :=  r[k]; 

TAB[ai,l]  :=  TAB[al,f]  -  TAB[af,al]  x  ref; 
end; 
end; 
end; 
end; 

b  :=  0;  al  :=  1; 

comment  The  following  segment  of  the  program  performs  a  column-by-column 
Gaussian  reduction  of  the  matrix  associated  with  the  reference 
equations,  forming  an  upper  and  a  lower  triangular  matrix  into 
the  array  P  .  (Each  diagonal  element  of  the  lower  triangular 
matrix  is  one. )  Interchanges  of  rows  take  place  so  that  the 
largest  pivot  in  each  column  is  employed.  It  is  assumed  that 
b-1  columns  have  already  been  decomposed.  If  the  matrix  is 
not  of  full  rank,  the  exit  insufficientrank  is  taken,  and  it 
is  left  up  to  the  user  to  determine  if  the  given  overdetermined 
system  can  be  solved  exactly. ; 
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body: 

tO  b;  11  :  =  b+1;  101  :=  b-1; 
for  i  :=  b  step  1  until  n  do 
begin 

t  :=  ix[!0] ; 

trisolv  (if  i=b  then  0  else  b, 1,101, 0>J-1, j, k, P[10, r[k] ] , 
if  r[j]«n  then  d[i]  else  A[l,r[j]],P[k,r[j]],l,ipr); 
trisolv  (10,l,n, 0,101, j,k,P[10,r [k] ], 

r[j]*n  then  d[l]  else  A[l,r[j] ],P[k,r[j]],l,ipr); 
ref  :«  0; 

for  j  :=  to  step  1  until  n  do 
begin 

t  :=  P[10,r[j]]j 

if  ref  <  abs  (t)  then  begin  ref  :«=  abs  (t);  s  :=  t;  k  :=  j;  end; 
end; 

if  ref  =  0  then  begin  j  :=  1;  _go  _to  singular;  end; 

10  =  n  then  po  to  decorapbranch[al] ; 
j  :=  r [k] ;  r[k]  :=  r[10];  r[10]  :=  j; 

for  j  :=  11  step  1  until  n  do  P[10,r[j]]  :=  P[10,r[j] ]/s; 

101  :*  10;  10  :=  11;  11  :=  11+1; 
end; 

singular: 

for  i  :=  0  step  1  until  n  do  refset[i]  :=  ix[i]; 

£0  to  failures[j]; 
return: 

comment  Solve  for  the  lambdas. ; 
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trisolv  (b,l,n,0, j-1, j,k,sv[k],  _if  r[j]-n  then  -1  else  0; 

P[k,r[j]],l,ipr); 

trisolv  (n,-l,0, j+l,n, j,k,lam[k],sv[j],P[k,r[j]],P[j,r[j3],ipr); 

comment  Compute  epsilon  for  the  reference  subsystem  of  equations. ; 
t  :=  0; 

for  i  :=  0  step  1  until  n  do  t  :=  t+abs(larafi]); 
eps  :=  l/t; 

comment  Each  new  value  of  eps  must  be  greater  than  the  previous  one. 

If  this  is  not  so,  the  solution  may  have  been  "overshot".; 
if  eps  <  lasteps  then  go  to  ed; 
lasteps  :=  eps; 

comment  Solve  for  the  vector  x,  the  Chebyshev  solution  of  the  reference 
subsystem  of  equations.; 

for  i  :=  0  step  1  until  n  do  xr[i]  :=  sign(lani(i ] )  X  eps; 
trisolv  (0,l,n,0,i-l,i,  j,w[j],xr[i],P[i,r[j]  ],P[i,r[i]],ipr); 
trisolv  (n,-l,0,i+l,n,i, j ,x[r ( j]  ] ,w[i  ] ,P[i,r  [  j ]  ]  ,l,ipr); 

comment  x[n]  should  be  -1  .  It  can  be  used  to  purify  eps  and  the  other 
components  of  x  . ; 
ref  :=  -x[n]; 

for  i  :=  0  step  1  until  nl  do  x[ i 3  :=  x[i]/ref; 
eps  :=  eps/ref; 

comment  For  each  index  ix[n+l],...,  ix[m-l]  compute  the  residual 
A[ix[j],o]  x  x[o]  +  ...  +  A[ix[j],n-l]  x  x[n-l]  -  d[ix[j]]  . 

If  the  largest  of  these  in  magnitude  is  not  greater  than  eps, 
go  to  itr  to  refine  the  vector  x,  for  it  may  be  the  Chebyshev 
solution  of  the  full  system. ; 
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resid  (ipr); 

if  ref  <  eps  then  go  to  itr; 
ovr: 

k  :■  ixfal]; 

comment  The  following  linear- system  solution  is  computed  in  order  to 
determine  which  equation  is  to  be  dropped  from  the  reference 
set  of  equations.  ; 

trisolv  (0,l,n,0,i-l,i, j ,w[ j 3 ,  r[i]=  n  then  d[k] 

else  A[k,r[i]],P[j,r[i]],l,ipr); 
trisolv  (n,-l,0,i+l,n,i, j,w[j],w[i],P[j,r[i]],P[i,r[i]],ipr); 
comment  s  is  the  sign  of  the  residual  with  greatest  magnitude.  Find 

the  largest  of  the  ratios  w[k]/lam[k]  X  s  .  If  any  component 
of  lam  is  zero,  the  exit  zerolambda  is  taken.  ; 
ref  :=  lam[n];  b  :=  n; 

if  ref  =  0  then  begin  j  :  =  2;  jjo  to  singular;  end; 

ref  :=  w[n]/ref  x  s 

for  j  :=  0  step  1  until  nl  do 

begin 

t  :=  lam{j]; 

if  t=0  then  begin  j  :=  2;  _go  _to  singular;  end; 
t  :=  w[j]/t  x  s; 

if  t  >  ref  then  begin  b  :=  j;  ref  :=  t;  end; 
end; 

comment  Form  a  new  reference  subsystem  by  exchanging  the  ix[al]-th 
and  ix[b]-th  equations. ; 

ix[ai]  :=  ix[b];  ix[b]  :=  k;  al  :**  1;  _go  ^  body; 

ed: 
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comment  Restore  the  previous  reference  subsystem. ; 
eps  :=  lasteps;  al  :=  2; 

j  :=  ix[al];  ix[ai]  :=  ix[b];  ix[b]  :=  j;  go  to  body; 


lasteps  :=  0;  cnt  :=  0; 
comment  Iteratively  refine  the  vector  x; 


cnt  :=  cnt  +  1;  if  cnt  >  10  then  go  to  insufficientrank; 

cnorm  :=  snorm  :  =  0; 

for  t  :=  0  step  1  until  n  do 

begin 

k  :=  ix[i]; 

t  :=  abs  (x[i]); 

if  snorm  <  t  then  snorm  :=  t; 

rv[i ]  :=  -ip2  (j,0,n,x[j],  _if  j*n  then  d[k]  else  A[k,j],  -xr[i]); 


trisolv  (0,l,n,0,i-l,i,j,rv[j)iPv[i],P[i,r[j]],P[iir[i]],ip2)i 
trisolv  (n,-l>0,i+l,n,i,j,w[r[j]3,rv[i],P[i,r[j]],l>ip2); 
for  i  :=  0  step  1  until  n  do 


begin 


s  :=  w[i ] ; 
x [ i ]  :=  x[i]  +  s; 
s  :=  abs  (s); 

if  cnorm  <  s  then  cnorm  :=  s; 


if  cnorm/snorm>  eta  then  go  to  ilp; 


u 


o  v-  >  .  • 


comment  eta  is  to  be  preset  with  a  small  positive  multiple  of  the  largest 


positive  single-precision  machine  number  u>  having  the  property 
that  1+to  =  1-u)  =1  in  single-precision  arithmetic.  The  small 
multiple  will  depend  upon  the  peculiarities  of  the  machine's 
rounding  process  and  will  have  to  be  empirically  determined. ; 
ref  :=  -x[n] 

for  i  :=  0  step  1  until  nl  do  x[i]  :=  x[i]/ref; 
eps  :=  eps/ref; 

comment  Determine  whether  a  Chebyshev  solution  has  been  found.  If  any 
residual  is  greater  in  magnitude  than  eps  while  eps  is  smaller 
than  a  value  produced  from  an  earlier  refinement,  give  up,  print 
a  warning,  and  return  the  best  x  computed  thus  far. ; 
resid  ( ip2 ) ; 

•  if  ref  <  eps  then  finished  :=  true 
else  if  eps  >  preveps  then  finished  :=  false 
else  begin  comment  Print  out  "DOUBTFUL  SOLUTION"; 

go  to  skip;  end; 

preveps  :=  eps;  refset[n]  :=  ix[n]; 
foi  i  :=  0  step  1  until  nl  do 
begin 

refset[i]  :=  ix[i]; 

h [ i ]  :=  x [i 3 ; 


end; 

if  -i  finished  then  go  to  ovr; 


skip: 

epz  :=  preveps; 
end  Chebyshev; 
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9*  Sample  Runs 


The  output  reproduced  in  the  appendix  was  produced  by  four  programs 
implementing  the  exchange  method.  At  each  exchange  step  the  reference 
set,  value  of  values  for  the  x^,  and  the  non-reference  residuals 

were  listed  followed  by  the  equations  to  be  switched  in  the  next  exchange. 
Upon  termination,  a  count  of  exchanges  and  solution  refinements  (where 
applicable)  was  printed  along  with  the  computation  time  required 
(print  time  excluded).  The  computed  Chebyshev  solution  for  the  full 
system  was  then  printed  followed  by  the  final  reference  set  and  a  list 
of  all  residuals. 

A  common  data  system,  Ax  ■  d,  was  given  to  the  four  programs. 

The  matrix  A  consisted  of  the  17  X  9  Hilbert  matrix  segment 

ai,j  =  i+I+I  16 5  8)  . 

The  right-hand  vector  d  had  components 

dj_  ■  i  (i=0, . . . ,  16)  . 

Output  A  was  produced  by  a  version  of  the  program  given  in 
section  8  using  double-precision  arithmetic. 

Output  B  was  produced  by  a  program  using  the  techniques  out¬ 
lined  in  section  4.  This  program,  however,  based  its  computation  on 
the  matrix 


at 

X1 

•  •  • 

at 

i 

n+1 

-d. 

X1 

•  •  • 
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rather  than  on  the  matrix  C  .  This  permits  the  initial 

_3  ? 

y  *  0(n2) 

operations  for  the  calculation  of  the  X^  to  be  saved ,  for  the  last 
column  of  B  satisfies 

/ 


Now,  however, 


must  be  computed  separately  at  each  exchange.  Note  that,  on  the  sample 
data,  this  program  has  failed  to  recognize  the  terminal  reference  set, 
giving  the  wrong  answer. 

The  suggestion  has  been  made  that  the  exchange  method  be  imple¬ 
mented  using  Jordan  elimination  techniques,  but  that  a  section  of  code 
be  provided  to  clean  up  the  solution  once  it  has  been  attained.  Output 
C  was  produced  by  such  a  program.  Clean-ups  were  carried  out  in  double- 
precision.  Since  this  program,  just  as  program  B,  failed  to  recognize 
the  final  reference  set  at  the  first  encounter,  the  clean-up  section 
was  called  upon  twice  for  the  given  data  set  -  once  to  put  the  program 
back  on  the  right  track,  and  once  for  the  final  solution  refinement. 

By  good  fortune  the  final  reference  set  was  recognized  the  second  time 
around. 
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Output  D  was  produced  by  a  B5500  Burroughs  Extended  Algol  version 
of  the  procedure  given  in  section  8. 

10.  Double-Exchange  Algorithm 

Instead  of  introducing  one  vector  into  the  reference  set,  we  con¬ 
sider  the  problem  of  introducing  two  vectors  simultaneously.  (What 
follows  can  easily  be  generalised  to  the  problem  of  introducing  several 
vectors  simultaneously. ) 

Without  loss  of  generality,  we  assume  that  A^, . . . ,  An+^  form  a 

reference  set.  Let  \  be  such  that 

1  n+i 

n+1 

£  VSt  =  0 

k=l  k  * 


under  the  normalization 


r.+l 

X.  v’k  *  -1  • 

k=l 

Then 

n+1 

€  =  1/  ]r  l\l  >  °» 

and  if  x  is  the  Chebyshev  solution  for  this  reference  subsystem, 
sgn(\k)  =  sgn(rk(x))  for  k=l,...,  n+1  . 
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For  ease  of  notation  we  write 


=  sgn(ri(x))Ai  for  all  i 
Tk  =  sgn(rk(x))\k  =  Ixkl  for  k=l,...,  n+1  . 

Thus 


I  tA  ■ 0 

k-i  K  K 


and 


n+1 

“  V  I  Tk 

k=l 


We  assume  that 


|r  (x) |  >  |r  (x)|  >  e 
*  1  a2 


for  some  a2  >  n+l  •  Since  Bn+1  have  rank  n 

exist 


9  •  •  •  } 
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’n+1 


and 


u(2) 
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so  that 


n+1 


k=l 


(J) 


B,  for 


3=1,  2  . 


The 


will  be  unique  if  we  also  demand  that 


n+1 


for  j=l,  2 
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We  wish  to  find  rows 


n+l})  to  exchange  with 


A  ,A  in  order  to  form  a  reference  set  with  a  greater  reference 
al  a2 

deviation  *'  .  Associated  with  this  will  be  a  reference  subsystem 
Chebyshev  solution  x'  .  Demanding  suitable  agreement  between  the  signs 
of  rk(x)  8111(1  (x ' ) ,  we  may  use  the  characterization  theorem  of 
section  2  to  determine  0^  and  02  .  Viz. ,  we  ask  for  numbers 
and  y2  such  that 


n+l 


Vl\  +  Y2Bff2  +  ^  (t1  *  -  0 


.(2)' 


with 


t!  =  y,  >  o  for  j-1,  2 
j  J 


=  (t< 


Vi 


(1) 


v2^2))  >  o 


for  i=l,...,  n+l 


and  for  two  indices 


The  normalizations  of  the  have  been  chosen  so  that 
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We  wish  to  choose  y^Yg  under  the  above  constraints  so  as  to  maximize 
c'  •  This  is  equivalent  to  determining  the  minimum  of 


Since 


n+1 

I  T, 

k=l 


is  fixed,  and  (as  can  easily  be  shown) 


n+1 


ej  3  L  ^ 

J  k=l 


E  JJ)  -  1  >  0  (j=l,  2), 


we  wish  to  determine  y^Yg  >  0  so  that 


Vi  +  Va 


is  maximized  subject  to 


+  V242)  <  %  for 


n+1 


35 


1  •' 


,  V  *■*. . 


>-r 


o  ?.<  \ 

h  v  r-  _  **  T  '  ' 

*  v  *s»^v  ^  . •  r_5  ^  •- *  * !  1  - 


This  is  a  standard  linear  programming  problem.  Note  that  the  single¬ 
exchange  algorithm  can  be  expressed  as  the  above  problem  with  the  addi¬ 
tional  constraint 


y2 


=  0  . 


Thus  the  e'  of  the  double -exchange  can  be  no  less  than  the  c'  given 
by  the  single  -exchange  of  section  3-  Note  further  that  conditions  1 
and  2  of  section  3  do  not  appear  in  the  development  of  the  double - 
exchange. 

Computation  can  be  simplified  by  considering  the  dual  to  the  above 

linear  programming  problem.  We  introduce  the  surplus  variables 

z  .o t  z  and  minimize 
n+2  n+3 


n+1 


T  2, 


i=l 


subject  to 


zi  >  0  for  all  i, 


n+1 


£  ^l)zr  ‘  Zn+2  c  *1' 


K— 1 


and 
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.(a). 


I  =  e, 

k=l  K 


k  "  "n+3  w2 
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If  either  surplus  variable  is  nonzero  in  the  solution,  then  B 

®] 

and  B  cannot  simultaneously  be  introduced  into  the  reference  set. 


The  correct  single-exchange,  however,  is  then  readily  obtainable  from 
the  dual  problem  solution. 

In  section  16  are  presented  some  timing  results  from  a  program 
implementing  this  algorithm.  Comparing  these  results  with  those  from 
the  single-exchange  implementations  of  sections  13-15,  we  see 
that  the  extra  effort  involved  is  not  paid  for  by  a  net  reduction  in 
time.  Also  we  have  observed  that  in  practice  rather  less  than  half 
of  the  exchange  steps  carried  out  permit  the  simultaneous  switching 


of  two  reference  equations. 


11.  Computational  Comparisons  of  Variations  for  the  Exchange 


In  the  procedure  given  in  section  8,  the  non-reference  equation 


chosen  to  enter  the  reference  system  at  each  exchange  was  the  a  , 
whose  residual  satisfied 


(a)  |r  (x)|  =  _  maX  \|r.(x)| 

'  'a  ijEl  reference  setl  i 

{  indices  J 


According  to  the  theory,  however,  the  exchange  method  will  converge  so 
long  as  the  reference  deviation  after  each  exchange  exceeds  the  refer¬ 


ence  deviation  before.  And  for  this  to  be  true,  it  is  sufficient  only 
that  a  satisfy  |r^(x)j  >  |e|  (conditions  1  and  2  given  section 
3  being  assumed  always  to  hold). 


I 


<4  V' 


K.«wr»  -jws?r  SKXP 


Alternate  versions  of  the  procedure  presented  in  section  8  were 
prepared  for  Stanford's  B5500  wherein  the  few  statements  determining 
a  according  to  (a)  were  changed  for  statements  implementing  other 
selection  rules.  The  unaltered  procedure  and  the  alternates,  together 
with  an  implementation  of  the  double-exchange  method  described  in 
section  1C,  were  run  on  random  systems  of  equations  of  several  sizes. 
Averages  of  times  required  and  number  of  exchanges  made  are  given.  Note 
that  the  procedure  of  section  3  gave  the  most  favorable  times. 

12.  The  Data 

Data  for  the  comparison  runs  was  generated  by  a  procedure  written 
in  Burroughs  Extended  Algol.  The  procedure  produced  a  matrix 

i=o,...,  m-1 

A  —  [fl.  . ]  i  -I 

ij  j— o, * . . ,  n-1 


and  a  vector 


d  *  [di]  i=o, . . . ,  m-1 

each  of  whose  elements  had  the  form  5X1),  where  5  was  a  pseudo¬ 
random  variable  distributed  approximately  uniformly  in  the  interval 
[o,  +1],  as  computed  by  the  mixed  congruential  method 

f=o 

*0 

)  ?n+l  =  (sl1  '  5)?n  +  2115271^9  mod  22? 
for  n  >  1  , 
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I 

r 

i 

i 

L 

l 

t 

L 

I 


L 

L 

U 


and  T|  was  chosen  pseudo-randomly  from  among  the  numbers 

+1,  +8"1,  +8’2,  +8’5,  -1,  -8'1,  -8“2,  -8"5  . 

Every  decision  rule  was  applied  to  ten  system,  each  of  m  equations 
in  n  unknowns,  where 

(m,n)c{  (10,4), (20,4), (30,4), (40,4),  (20,9), (30,9), (40,9), (30,19)}  • 

13*  Selection  of  the  Equation  with  Largest  Residual  Magnitude 

The  procedure  given  in  section  8  produced  the  following  statistics 
e  mean;  a  =  standard  deviation): 


Time  Required  (Seconds) 


nN^n 

4 

9 

19 

10 

H=0.677 

O=0.110 

20 

H=1.079 

^=4.043 

0=0. 142 

o=0. 850 

30 

1**1.246 

^=5.947 

^=28.620 

o=o. 236 

o=l. 170 

0=6.802 

40 

HFl. 558 

^=7 . 265 

o=0. 266 

0=1-740 
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Number  of  Exchanges 


"Yj 


10 

H=3-40 

a=l. 56 

20 

H=5-90 

H=9-10 

o=l. 8l 

o=3-05 

30 

^=5-90 

n= 13-40 

H=l6.80 

0=2.21 

o=3-64 

0=5-21 

40 

^=6.70 

H=l4.60 

0=2. 10 

o=5- 16 

14.  Selection  of  the  First  Suitable  Equation  Found 


The  first  variant  program  examined  each  non-reference  equation  in 
turn  until  one  was  found  whose  residual  magnitude  exceeded  the  reference 
deviation.  That  equation  was  selected  for  introduction  into  the  ref¬ 
erence  system.  Statistics  for  this  variant  follow. 
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15*  Selection  So  As  to  Give  the  Greatest  Reference  Deviation  Increase 


Given  any  non-reference  row 


A  for  which 
a 


I rcy ( x ) I  > 


solve 


Then,  if  0  is  such  that 


sgn (e)  sgn(ra(x))^ 

X3 


is  maximal,  A  must  replace  A.  in  the  reference  set.  The  new 

“  e 

\'s  can  be  computed  as  follows: 


S  =  Xp/*p 


=  *i  -  -  Xp  (*#> 

H3 


Then 


'I  |r#(x)|  *  (1  -i^l)  lei  , 


where 


K  = 


n+1 

k 


ho 


Using  these  results,  a  variant  of  the  procedure  given  in  section  8 
was  prepared  in  which  the  non-reference  equation  selected  to  enter  the 
reference  system  at  each  exchange  was  that  one  which  would  give  the 
greatest  value  to  | e  *  I 


Time  Required  (seconds) 


m\n^ 

4 

9 

19 

10 

4=0.821 

0=0. 187 

20 

4=1-315 

0=0.327 

4=5-900 

o=l-886 

30 

4=1. 528 

4=9-798 

4=42.481 

o=0.313 

0=2.423 

0=7-921 

40 

4=2. 134 

o=0. 465 

4=l4. 685 

0=3-825 

4l 


*• 

J.  if*. it 


* 


V' 


& 
f v 


/  . 

.  v  » 

V 


Number  of  Exchanges 


m\^r 

1  4 

9 

19 

10 

H=3.50 

o=l*  56 

20 

l*=4. 60 

(*=8.10 

0=1*7*+ 

o=2. 84 

50 

^=5.90 

(*= 10.40 

(*=lo. 20 

o=l. 04 

o=2.20 

0=3.03 

40 

(*=5« 10 

(*= 15.40 

0=1.70 

p=5.23 

lb.  Double -Exchange  Algorithm 

Time  Required  (seconds) 


_ 9  19 


10 

(*=0-900 

0=0.147 

20 

(*-1.258 

o=0.215 

n=4.557 

0=0.836 

30 

(i=1.442 

0=0.271 

(*=6. 487 

o=0. 951 

(*=36.650 

0=8.179 

40 


1.912  n=9.4lj 
0=0.677  0=1-507 


(An  exchange  cycle  consisted  of  the  simultaneous  switching  of  two 
equations  where  possible.  Otherwise  it  consisted  of  a  standard  single¬ 
exchange  . ) 
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****  HILBERT  DATA  ***** 

17  EQUATIONS  IN  9  UNKNOWNS 


EXCHANGE  ALGORITHM  IN  DOUBLE-PRECISION 


REFERENCE  SeTi 

0  2  11  1516  3  849 

1,65566# 11074,00287,19989,63#  -3  ■  EPS 

3,97047# 60096, 43108,06759, 279  3  ■  X[  0] 

-2.70355,98109,22439,41451,339  5  ■  X[  1J 

4.55974,74511,68592,81444,739  6  a  XC  23 

-3.26294,80229,02832,31393,358  7  ■  X[  33 

1.20427,47950,72163,13981,29#  8  -  X[  43 

-2. 48030# 00464, 27790# 08801 #459  8  a  Xf  53 

2. 87782#84818#63886, 24654*619  8  ■  Xf  63 

•1.75780#5371O#6389O#56612#9O8  8  ■  XC  73 

4.39490,84805# 13472*43622*049  7  ■  X[  83 

•2.57761,12267,27838,37828,099  -3  ■  RESIDUALt  103 
l.C7571#3597B#74505#82246#038  -3  a  RESIDUALt  63 
1.00245,84291# 63795,29485,838  -2  ■  RESIDUALt  123 
1.95676,92968,81866,18250,259  -2  «  RESIDUALt  133 
2.48292,91781,18352,39794,149  -2  a  RESIDUALt  143 
1.98447,14020,05852,36957,669  -2  a  RESIDUALt  153 
3.46169,10322,94211,82277,299  -3  •  RESIDUALt  73 
EXCHANGE  EQUATION  11  WITH  EQUATION  14 

REFERENCE  SeTi 

0  2  14  1  5  16  3  8  4  9 

3.29205,02104,33056,37236,26#  -3  a  EPS 

5.40554,73215,77523,85646,17#  3  ■  X[  03 

•3.59146,43813,84385,18429,92#  5  ■  X[  13 

5.93749,90812,11901,92471,069  6  a  Xt  23 

-4,17859,65144,44931,34071,229  7  a  XC  33 

1.52049,75097,12577,52825*279  8  a  XC  43 

-3.09350,90379,71150,66375,699  8  a  X[  53 

3.55123,22422,42775,98985,69#  8  a  Xt  63 

-2.14895,91603,99356,68537,189  8  a  XC  73 

5.32811,40323,12399,58580,08#  7  a  XC  83 

-1.04956,72051,91559,43804,889  -2  ■  RESIDUALt  103 
-2.27328*65222,76318,77735,809  -4  a  RESIDUALt  63 
•1.08864,60005,98963,99208,988  -2  a  RESIDUALt  123 
-3.91191,80148,65456,42774,218  -3  a  RESIDUALt  133 
-1.34090,31468,42742,36611,82#  -2  a  RESIDUALt  113 
5.51335,12173,87940,65038,29#  -3  a  RESIDUALt  153 
4.64657,18263,68794,95518,65#  -3  a  RESIDUALt  73 
EXCHANGE  EQUATION  9  WITH  EQUATION  11 

REFERENCE  SET» 

0  2  14  1  5  16  3  8  4  11 

5.30006,47585,98979,14408,70#  -3  a  EPS 

6.27879,92051,09165,06026,78#  3  «  XC  03 

i 


-4.C9612,  36199,45232,02513,46#  5  a  Xt  i) 

6.67733,07627,41189,48351,34#  6  ■  Xt  23 

-4.64706,76438,19307,65657,34#  7  ■  X[  3] 

1 ,67565,40534,98819,77968,69#  8  ■  Xf  4] 

-3.38355,28056,64260,96258,72#  8  a  Xt  53 

3.85958,25436,98984,77619,29#  8  *  Xt  63 

-2.32292,40280,66317,48340,80#  8  ■  Xt  73 

5.73258,  79224,06205,97422*59#  7  ■  Xt  83 

-2.61547,83860,89942,20876,97#  -3  ■  RESIDUAL!  103 
-4.40861,93560,47416,28225,15#  -3  a  RESIDUAL!  63 
-4.10282,78390,15665,83395,49#  -3  ■  RESIDUAL!  123 
6.41483,09155,73710,43478,49#  -4  ■  RESIDUAL!  133 
2, 48599, 96575, 45978, 62499, 50#  -3  ■  RESIDUAL!  93 
5.20435,98784,03175,26575,97#  -3  ■  RESIDUAL!  153 
2.53687,03218,72219,94661,28#  -3  a  RESIDUAL!  73 

TERMINATION 

NUMBER  OF  EXCHANGES  MADE  WAS  2 
TIME  IN  SECONCi  ■  5,48 

SOLUTION  vectori 

6.27879,92051,09165,06026,78#  3  a  Xt  03 

-4. 09612,36199,45232,02513,46#  5  a  Xf  13 

6.67733,07627,41189,48351,34#  6  a  Xf  23 

-4.64706,76438,19307,65657,34#  7  a  Xt  33 

1.67565,40534,98819,77968,69#  8  a  Xt  43 

•3.38355,28056,64260,96258,72#  8  a  Xt  53 

3.85958,25436,98984,77619,29#  8  a  Xt  63 

•2.32292,40280,66317,48340,80#  8  a  Xt  73 

5.73258,79224,06205,97422,59#  7  a  Xt  83 


REFERENCE  SET  I 

C  2  14  1  5  16  3 

RESIDUALS  t 

5. 30006# 47585# 99232# 206 39# 30# 
-5. 30006# 47585# 98899# 13948# 56# 
5,30006#  47585*99232#  206  39#  30# 
-5, 30006# 47585# 98677# 09488# 07# 
5, 30006# 47585# 99232# 20639# 30# 
*5, 30006# 47585# 991 76# 69524# 17# 
-4. 40861, 93560# 474 16, 28225# 15# 
2, 5  3687, 03218, 7221 9# 94661 #28# 
5, 30006# 47585# 96954# 6506 3# 68# 
2. 48599# 96575# 45978# 62499# 50# 
-2. 6 1547# 83860# 89942# 20876# 97# 
•5, 3 0006# 47585# 99232# 206 39# 30# 
-4. 10262..  78390#  15665,83395,49# 
6 ,4 148 3# 091 55# 737 10,43478 #49# 
5,30006,47585,98899# 13948,56# 
5, 204 35# 98784# 03175# 26575# 97# 
-5, 30006# 47585# 99287# 71754# 42# 


8  4  11 


•3 

■ 

RESIDUAL! 

0) 

•3 

■ 

RESIOUALt 

M 

•3 

8 

RESIDUAL! 

21 

•3 

■ 

RESIDUAL! 

3) 

•3 

■ 

RESIDUAL! 

41 

•3 

u 

RESIDUAL! 

51 

•3 

w 

RESIDUAL! 

61 

■3 

u 

RESIDUAL! 

71 

•3 

w 

RESIDUAL! 

81 

•3 

a 

RESIDUAL! 

91 

■3 

■ 

RESIDUAL! 

101 

3 

a 

RESIDUAL! 

111 

3 

a 

RESIDUAL! 

121 

4 

■ 

RESIDUAL! 

131 

3 

a 

RESIDUAL! 

141 

3 

a 

RESIDUAL! 

151 

3 

a 

RESIDUAL! 

16] 

s 


>«**m**v 


*****  HILBERT  DATA  *****  g 

17  EOUAT IONS  IN  9  UNKNOWNS 
tableau-jordan  ALGORITHM 
COMPUTATION! 


REFERENCE  SeT! 

0  2  11  1  '5  16  3  8  4  9 

EPS  ■  1 ,65552 3068334*0 3 

XI  0]  ■  3,967545797508*03 

X[  1)  ■  “2,701867054788*05 
XC  21  ■  4,557296022678*06 

XC  33  ■  “3.261422535888*07 
X[  4]  «  1  ,203779412998*08 

XC  5]  ■  -2,479392868428*08 
XC  6]  *  2,876883648618*08 

XC  7]  ■  -1.757363248278*08 
XC  83  ■  4,393718637028*07 

RESIOUALC  103  «  “2,593994140634-03 
RESIDUALC  63  -  1.129150390638-03 

RESIOUALC  123  ■  1.001358032234-02 

RESIDUALC  133  *  1 .957321 166994-02 

RESIOUALC  143  *  2*487564086918-02 

RESIDUALC  151  ■  1.987075805664-02 

RESIDUALC  73  *  3.509521484384-03 

EXCHANGE  EQUATION  14  WITH  EQUATION  ll 


REFERENCE  SET! 

0  2  14  1  5  16  384  9 

EPS  a  3.283101470788-03 
XC  03  «  5,394102722128*03 

XC  13  *  -3,584534690988*05 
XC  23  *  5.926966466308*06 

XC  33  «  “4.171726764028*07 
XC  43  «  1,518164465358*08 

XC  53  ■  “3.089051809824*08 
XC  63  *  3.546402234558*08 

XC  73  ■  “2.146187334578*08 
XC  83  ■  5.321569977128*07 

RESIDUALC  103  a  -1.046752929698-02 
RESIOUALC  63  a  2.134230468754-04 
RESIDUALC  123  ■  -1 .065063476568-02 
RESIDUALC  133  «  -3.570556640638-03 
RESIOUALC  111  a  -1.312255859384-02 
RESIOUALC  153  a  5.706787109408-03 
RESIDUALC  73  a  4.852294921888-03 
EXCHANGE  EQUATION  11  WITH  EQUATION  9 


REFERENCE  SeTi 
0  2  14 


16 


8 


4  11 


1 


! 


r 

I 

I 


EPS 

S 

5.274902555224-03 

X! 

0] 

=  6.259571972104403 

X! 

1] 

■  -4.004668078720*05 

X! 

2] 

»  6,660130004704*06 

X! 

3] 

■  -4, 635957755864*07 

X! 

4] 

*  1.671911561594*08 

X! 

5  ] 

■  -3.376452397544*00 

X! 

6] 

a  3, 851934333774*08 

X! 

7] 

*  -2.318558082030*00 

X! 

8] 

■  5.722328674504+07 

RESIDUAL!  10]  s  -2.471923428130-03 
RESIOUALC  6]  ■  -4.028320312500-03 
RESIDUAL!  12]  »  -4,150390625004-03 
RESIDUAL!  13]  ■  A , 544921 875004-04 
RESIDUAL!  9)  =  2 . 7 1 606445 3 1 3P-0 3 

RESIDUAL!  15)  ■  5 . 40 1 6 1 1 323 1 34-0 3 

RES iOUAL I  7]  «  2.868652343754-03 

EXCHANGE  EQUATION  15  WITH  EQUATION  14 


REFERENCE  SeTi 


0 

2  15  1  5 

EPS 

■ 

5.286956166464-03 

X! 

0] 

«  6.266837272300*03 

X! 

11 

«  -4,089030849654*05 

X! 

2] 

a  6,666739045000+06 

X! 

3] 

■  -4,640263866400+07 

X! 

4] 

■  1.673374064974*08 

XI 

5] 

a  -3,379248217180*08 

XI 

6] 

*  3.854966854230*08 

XI 

7] 

«  -2.320300230044*08 

XI 

8] 

«  5.726446347104*07 

RESIDUAL!  101  « 
RESIDUAL!  6]  « 
RESIDUAL!  121  ■ 
RESIDUAL!  13]  * 
RESIDUAL!  9]  « 
RESIDUAL!  14]  « 
RESIDUAL!  7]  . 


-2.716064453134-03 
-4 ,4250404281 34-03 
-4,089355468754-03 
4,577636718754-0' 
2. 30O37lO9375R'  3 

5.249023437*04'-  )3 
2, 807617147508-03 


3  6 


3 


50 


TERMINATION 

NUMBER  OF  EXCHANGES  WAS 


; 


I 

l 

|  time  in  seconds  =  ?.ifl 


SOLlTiriN  v E C TOR  I 


r 

xt 

01 

=  6.26683727?30P+03 

1 

xt 

1 1 

c  -4.0fl903084965P+05 

X[ 

23 

s  6.6AA73904500»*0A 

xt 

3] 

«  •4,A/j02A386A«0*+07 

xt 

4] 

*  1 ,67l  37<i0<)497M  +  0fl 

xt 

bl 

*  •3,370,wiH21  7  1UP+0H 

xt 

6] 

*  34flb«9668b42.3<**08 

I 

xt 

71 

«  “2. 3203002 J004H+08 

1- 

xt 

81 

a  !>.7?644  63«710fr+07 

r 


REFERENCE  SET l 
0  2  15  1 


5  16  3  8  4  11 


RESIDUALS  I 

Rt  0]  »  5.030549329638-03 
Rt  13  *  "5.422996449088-03 
Rt  23  ■  5.213819480858-03 
Rt  33  «  -5.371528525158-03 
Rt  43  *  5.212558797788-03 
Rt  53  *  -5.294128506538-03 
Rt  63  *  -4.351757061368-03 
Rt  73  *  2.577717493358-03 
Rt  83  *  5.295551415858-03 
Rt  93  *  2.4/11709321218-03 
Rt  10]  *  -2,676972399328-03 
RC  113  *  "5.355667924338-03 
Rt  123  *  -4.135860841088-03 
RC  133  *  6.373025899608-04 
RC  143  *  5.322239438788-03 
RC  15]  ■  5.242799588258-03 
RC  163  *  -5,260471501408-03 
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*****  hilaeht  data  ***** 

17  EQUATIONS  IN  9  UNKNOWNS 

tableau-joroan  algorithm  with  clean-ups 

COMPUTATION! 


REFERENCE  SfT  t 


0 

?  n  i  s 

EPS 

s 

1 

.65552306A33P-03 

X! 

0) 

8 

3. 9*7545797500*03 

X! 

1) 

a 

"2.7018670547 Du *05 

X! 

2) 

a 

4, 55729602267*'.  *06 

X! 

3) 

s 

-3,261422535886+07 

XI 

43 

a 

1 .2O377C4l29Vis*0R 

X! 

5) 

a 

-2.47939286Hu2r*0A 

X! 

6) 

a 

2.67688384AM  >UOR 

X! 

7) 

a 

"1.757363248279+08 

X! 

A] 

a 

4 , 3  9 .3  7 1  8  6  3  7  0 +  0  7 

lft 


RESIDUAL!  103  e  -  ? .  59  3994  1  406  3*»-0  3 
RESintlALC  63  s  1  .  1  2<>1  !iO  3906  3«»-0  3 
RESIDUAL!  123  *  1 . 00  \  3580322  3*-0? 
RESIDUAL!  13?  =  1 . 957 321 1 6690»-0? 

RESIDUAL!  143  *  ?.4875„.iCf6'V1 

RESIDUAL!  153  *  1 . 9*  7 075*0566'  -0? 
RESIDUAL!  73  *  3,S095?l 404 3BP-03 

EXCHANGE  FQUATION  1ft  WfTH  FQIJATION 


II 


4  9 


REFERENCE  SET: 


0 

2  14  1  5 

EPS 

s 

3 

.283101  4707«ia-03 

X! 

03 

8 

5. 3941027221 +03 

XI 

1 3 

a 

"3.58453469098^+05 

XI 

23 

a 

5. 9269664  6  6 30H +06 

X! 

33 

a 

"ft  .  1  7172676402^+07 

XI 

4] 

a 

1 .5181  *4465358+08 

XI 

53 

a 

-3.08O05t689P2H+0P 

XI 

6) 

a 

3,sl4640223ft65«  +  OA 

XI 

73 

a 

-2.1461P733457P+0A 

XI 

83 

a 

5. 3215699771  28  +  07 

RESIOUALC  103  «  -1.046752929698-02 
RESIDUAL!  63  *  ?  .  1 3623046B75P-04 

RESIDUAL!  1?3  »  -1 ,065Pft34765ftP-02 
RESIDUAL!  1  33  *  -  3 . 57O55664063P-03 
RESIDUAL!  11)  ■  -1.3122550593AP-O2 
RESIOUAL!  15)  ■  5.70678710940P-03 

RESIDUAL!  7)  *  4 , 65279492 1 88P-0 3 

EXCHANGE  EQUATION  ll  WITH  EQUATION  9 

reference  set i 

0  214  1  516  3  B  4  11 


/ 


Sr 


-WO 


EPS 

S 

5.274902555228-03 

X! 

0] 

■  6.25957!97210W*03 

Xt 

1  ] 

s  -4, 08466807872 8 *05 

X! 

2] 

*  6.66013000470*406 

X! 

31 

«  -4, 6 3595 7 75586ft 407 

X! 

4] 

■  1.67191156159*408 

X! 

5] 

a  -3.37645239754*408 

X! 

6] 

«  3.85193433377*408 

X! 

7] 

*  •?«  31855808 20 3*408 

X! 

8] 

■  5.72232867450*407 

RESIDUAL!  10]  s  -2.4719238281 38-03 
RESIDUAL!  61  e  -4.020320312508-03 
RESIDUAL!  1?J  *  -A . 150390625008-03 
RESlDUALt  133  *  8 . 5*492 1 875008-0* 

RESIDUAL!  91  =  ? . 7 1 60ft A  45 3 1 38-03 

RESIDUAL!  15]  =  5 . A0 1  ft \ \  328 1 38-0 1 

RESIDUAL!  7]  »  ?. 8*8652 3A 3758-03 

EXCHANGE  EQUATION  15  WITH  EQUATION  1« 


REFERENCE  SETi 


0 

2  15  1  5 

EPS 

s 

5 

.286956166468-03 

XI 

0] 

2 

6.26683727230*403 

XI 

1  ] 

2 

-4.08903084965*405 

X! 

21 

e 

6.666739045008+06 

XI 

3] 

s 

-4,6*0253066*0*407 

XI 

4] 

= 

1 .67337*06*97*408 

XI 

51 

s 

-3.37924821718*408 

XI 

6] 

2 

3.85496605423*408 

XI 

7] 

= 

-2.32030023004*408 

XI 

8] 

= 

5.726446347108407 

RESIDUAL!  103  *  -2.716064453138-03 
RESIDUAL!  ft]  «  -A. 425048820138-03 
RESIDUAL!  l?3  »  -  A .089355460758-03 
RESIDUAL!  13]  =  4 .577636710758-0* 
RESIDUAL!  9]  s  2.380371093758-03 
RESIOUAL!  1*1  =  5.2A9021A375D8-03 
RESIDUAL!  73  =  2.007617187508-03 


dduple-hRecisidn  improvement 


reference  sf t * 

o  2  15  1  5  16  3 

5. 278  A*# ftOft9 3, 620 39, 81 16  3* 158 
6.26576,  18*67, 809*7, 90481* 138 
-4.00829, Afl« 47, 3608 A ,00  166,  *18 
6, 6 65 A  7,  1 A 70 8, 78 32 3 *27085 ,9 78 
-A, 6393**08 (J3P.7PSV0, 99856,  31* 
1 .67302.98250,58*17,63321,99# 
-3. 3 78b 3, 62 34 1,23599, 235 A  A, 098 
3,85414*1  3 358 *25 A  37,6400t»,688 
-2.31979,81579.50857.74466,808 
5,72519,98266,81440,61070,738 
-2.62825,04150,30675, 1 3559,698 


8  4U 


-3 

9 

EPS 

3 

m 

xt 

0) 

5 

n 

xt 

1] 

6 

w 

xt 

21 

7 

a 

xt 

3) 

8 

a 

xt 

4] 

6 

a 

xt 

5] 

8 

a 

Xt 

6] 

8 

a 

xt 

7) 

7 

a 

Xt 

61 

-3 

a 

RESlDUALt  10] 

51* 


5 nvwmn 


\ 


-4, 37417, 88687,2475  A.  19546#  1  3»  -1  -«  RrSlDUALt  61 

•4,04290,63798,48812. 16847, 2  t8  -3  *  RESIDUAL!  12) 
7,291  10,  *49988,26042,87869, 16*  -4  ■  RESIOUALl  133 
2,45568,04766, 19148,67774, IS*  -J  ■  RESIDUAL!  93 
5, 39490, 96191, 16761, 145M*0S«  -3  ■  RESIDUAL!  143 
2. 54498, 84901, 47729, 98658*  728  -3  >  RESIDUAL!  71 
EXCHANGING  EQUATION  15  WITH  EQUATION  14 

IMPROVEMENT  COMPLETE  -  SOLUTION  NOT  YET  ATTAINED 


REFERENCE  $ET| 

16  0  2  8  1  a  11  3  14  3 

EPS  «  5.335096772159-03 


X! 

03 

a 

6, 271517680108*03 

X! 

13 

a 

-4,091349270098*05 

X! 

23 

a 

6.660674302909*06 

X! 

33 

a 

-4, 641881114708*07 

XI 

43 

a 

1 .673847704968*08 

X! 

53 

a 

-3.380017444558*08 

X! 

63 

a 

3.855686242748*08 

XI 

73 

a 

-2.170657383018*08 

X! 

83 

a 

5.727176281408*07 

RESIDUAL!  103  ■  99072?65625*-03 

RESIDUAL !  63  ■  -4.516601562508-03 

RESIDUAL!  123  ■  -4.425048878139-03 
RESIDUAL!  133  ■  4 . 272460937509-04 

RESIOUALC  93  ■  2.197265625009-03 

RESIDUAL!  153  «  5.035400390639-03 

RESIDUAL!  73  «  7.197265625009-03 

double-precision  IMPROVF.MENT 
REFERENCE  SETj 

16  0  2  8  1  4  1  1  3  14  5 

5.30006,47585,99124,11354,729  -3  ■  EPS 

6,27879,92051,09148,09028,689  3  ■  X!  03 

‘-4.09612,36199,45217,87480,489  5  ■  X[  1] 

6.67733,07627,41163,91494,089  6  a  XI  23 

-4.04706,76438,19289,03505, 149  7  ■  X!  33 

1.67565,40534,98812,94935,719  8  *  XI  43 

-3.38355,28056,64247, 14797,909  8  ■  XI  53 

3,85958,25436,98969,13231,219  8  ■  X!  63 

-2.32292,40280,66308,18546,289  8  ■  XI  73 

5.73258,79224,06183,38542,499  7  ■  X!  83 

-2. 61547, 83860, 90108, 74792*34*  -3  ■  RESIDUAL!  103 
-4, 40861, 93560, 46417, 939  -3  ■  RESIDUAL!  63 
-4.10282,78390,15554,81165,749  -3  ■  RESIOUAL!  123 
6, 41483,09155,76685,99^4, 65*  .4  ■  RESIDUAL!  133 
2.48599,96575,45923,1 1  484  , 389  -3  ■  RESIDUAL!  93 
5. 20435, 98784,03230.77691, 099  -3  «  RESIDUAL!  153 
2.53687,03218,72775,05812,519  -3  ■  RESIDUAL!  73 

TERMINATION 


NUMBER  OF  EXCHANGFS  WAS  3 

NUMBER  OF  SOLUTION  RFFlNfMFNTS  WAS  2 

I 

- -  "  *  . . . .  . 


»>-■*  ■ms**'  j 


TIME  IN  SECONDS  «  7.?? 


SOLUTION  VECTOR! 


xc 

01 

■ 

6.278799205104*03 

xc 

M 

■ 

-4.096123619954*05 

X[ 

21 

■ 

6,677330762804*06 

xc 

31 

■ 

-4.64706764.3824  +  07 

xc 

41 

■ 

1.67*654053509*08 

xc 

51 

■ 

•3,383552805679+08 

xc 

61 

■ 

3. 859582543694*08 

xc 

71 

■ 

"2. 3? ?9 24028064*08 

xc 

81 

■ 

5,732587922404*07 

RFFCRFNCf!  S>  T  i 

18  0  ?  8  1  4  11 


RfStOUALSt 

Rt 

01 

■ 

5. 140681  667 14*-')  3 

Rt 

1) 

s 

-5.4  J(S94I06409«-01 

R  ( 

21 

■ 

5,180081  l744R«-03 

Rt 

31 

■ 

*5, 4069 330 76999-0 3 

R  [ 

41 

a 

5, 2037 2 05247? ^ -03 

Rt 

bi 

a 

-5, 3877^  1**??5  3  >-03 

Rt 

8) 

a 

-4,4891.3723409  .-03 

Rt 

71 

u 

?,  46245600966^-03 

Rt 

81 

9 

5. 230090727  302-0? 

Rt 

91 

9 

?.4pi  3742521  ?4,-03 

Rt 

101 

9 

“7.07011 79M70Hh-03 

Rt 

111 

9 

“5, 3s7tP?75?65'a“03 

Rt 

121 

9 

-4,15681  1  61849'--93 

Rt 

131 

8 

5, 90306844980, *-04 

Rt 

141 

a 

5,251 41 799559^-03 

Rt 

151 

a 

5. 1*5800396900(^03 

Rt 

1 8  1 

9 

“5. 3443361 48?pp-03 

*****  HILBERT 
17  EQUATIONS 
Q0LU8-RARTELS 


data  ***** 

IN  9  UNKNOWNS 
PROCEDURE 


D 


COMPUTATION! 


REFERENCE 
C  2 

EPS  ■  1.6 

XC  0) 

XC  1] 

XI  21 
XI  31 
XI  41 
XI  51 
XI  6] 

XI  71 
XI  8] 
RESIDUAL! 
RESIDUAL! 
RESIDUAL! 
RESIDUAL! 
RESIDUAL! 
RESIDUAL! 
RESIDUAL! 
EXCHANGING 


SC T  I 

11  1  s  16  3 

38740968008-03 
3.954765136728*03 
2.694271884778*05 
4.545869589848*06 
3.254013444538*07 
1.201271788238*08 
2.474611566238*08 
2.871708030238*08 
1.754394392028*08 
4.386710147658*07 
101 
61 
121 
133 
141  i 
151 
71  . 


-2.655029296888-03 

1.129150390638-03 

1.007461547858-02 

1.968765258798-02 

2.496337890638-02 

1.995086669928-02 

3.631591796888-03 


EQUATiON  11  WITH  EQUATION 


9 


14 


REFERENCE  SET: 

02  14  15  16  3049 

EPS  *  3,279427677398-03 

X!  0]  «  5,394266113288*03 

XI  11  *  -3.584942485358*05 
XI  21  =  5.927948390608*06 

X!  31  *  -4. 172558304108*07 
X!  41  =  1.518500985148*08 

X!  51  *  -3,089782861778*08 
X:  61  »  3.547276950468+08 

XI  7]  »  -2.146730646908+08 
XI  81  a  5.322939958218+07 
RESIDUAL!  101  a  - 1 , 046752929698-02 
RESIDUAL!  61  «  -6.103515625008-05 
RESIDUAL!  121  a  "1,089477539068-02 
RESIDUAL!  131  «  -3.875732421888-03 
RESIDUAL!  Ill  a  -1.324462890638-02 
RESIDUAL!  151  »  5 . 706 787 109408-03 

RESIDUAL!  71  a  4,943847656258-03 
EXCHANGING  EQUATION  9  WITH  EQUATION  11 


REFERENCE  SET! 


0 

2 

14  1  5  16 

EPS  a 

5. 

30913757  '208-03 

X! 

01 

= 

6, 27637304 ft 908+03 

XC 

11 

= 

-4 .094718295898+05 

4  11 


58 


X! 

23 

■ 

6.675261306609+06 

X! 

33 

■ 

*4,645779785749+07 

X! 

43 

■ 

1.675231013139+06 

X! 

53 

■ 

-3,382768370599+08 

X! 

63 

■ 

3.858754992099+08 

X! 

73 

a 

-2.322460488979+06 

X! 

83 

a 

5,731517357209+07 

RESIDUAL C  10]  ■  -2.593994140639-03 
RESIDUAL!  6]  ■  -4.394531250009-03 
RESIDUAL!  123  >  -4.150390625009-03 
RESIDUAL!  133  ■  6.406691406309-04 

RESIDUAL!  9]  ■  2.532956964369-03 

RESIDUAL!  153  ■  5.035400390639-03 

RESIOUAL!  7J  a  2.288816359389-03 

ITERATIVE  IMPROVEMENT 

REFINED  VALUES! 

EPS  ■  5.300064756599-03 

X!  03  >  6, 278799205 109+0 3 

XI*  13  ■  -4.096123619959+05 
XI  23  ■  6,677330762809+06 

XI  33  ■  -4,647067643829+07 
X!  4]  ■  1,675654053509+06 

XI  53  *  -3.383552805679+06 
XI  6]  «  3.859582543699+08 

X!  73  ■  -2.322924026069+08 
XI  83  ■  5.732587922409+07 

RESIDUAL!  103  ■  -2.676117987889-03 
RESIDUAL!  6]  «  -4.489137234699-03 
RESIDUAL!  123  ■  -4 . 15681 1618499-03 
RESIDUAL!  133  «  5.903068449809-04 

RESIDUAL!  9]  ■  2.421374252129-03 

RESIDUAL!  153  ■  5.158003969809-03 

RESIDUAL!  7]  «  2.462456009669-03 

TERMINATION 

NUMBER  OF  EXCHANGES  MAOfc  HAS  2 
NUMBER  OF  SOLUTION  REFINEMENTS  WAS  l 


59 


TIME  IN  SEC0N03  ■  3*10 

SOLUTION  VECTORI 
XI  03  ■  6. 278799 20510# 403 

X[  ll  ■  "ft  ,09612361995#405 
XI  23  ■  6,67733076280#406 

XI  33  ■  "ft  ,6ft70676ft382#407 
XI  ft 3  ■  l«67565ft05350#408 

XI  53  ■  -3.38355280567#40fl 
XI  63  ■  3.8595825ft3A9P*0A 

XI  73  ■  "2.32292ft02806M+08 
XI  83  ■  5.732587922ft0M40 7 


4  11 


t 


REFERENCE  SETl 


0 

2 

14 

RESIDUALS! 

RC 

0) 

■ 

5. 

RC 

n 

a 

-5. 

Rt 

23 

a 

5. 

Rt 

3] 

a 

-5. 

R[ 

4  3 

a 

5. 

Rt 

53 

a 

-5. 

RC 

63 

a 

-4, 

Rt 

73 

a 

2. 

Rt 

83 

B 

5. 

Rt 

93 

8 

2. 

Rt 

103 

8 

-2. 

Rt 

113 

8 

-5. 

Rt 

123 

B 

-4, 

Rt 

133 

8 

5. 

Rt 

14] 

a 

5. 

RC 

153 

a 

5. 

Rt 

163 

a 

•5. 

1  5  16 


8 


•4 


t 

> 


1 


l 

L 

[ 

L 


6l 


Unclassified 


Security  Classification  _ 


DOCUMENT  CONTROL  DATA  •  RAD 

ftMurttr  ilNiUlMllM  •/  Mite,  Ml  M  ateiMii  W  Mtital  mntlthm  mutt  M  M»M  «A—  te  wwll  »*A*W  te  iMMSM 


I.  OAIOINATIN  e  ACTIVITY  fC«*««te  «i«AM  «•  RBRRRT  •BCURITV  «  b  A44IPIS  A  TIM* 

Computer  Science  Department  Unclassified 

Stanford  University  ■>  snaur  " 

Stanford,  California  94305  _ --- _ 


S.  MI  MOAT  TITLB 

COMPUTATIONAL  CONSIDERATIONS  REGARDING  THE  CALCULATION  OF  CHEBYSHEV  SOLUTIONS  FOR 
OVERDETERMINED  LINEAR  EQUATION  SYSTEMS  BY  THE  EXCHANGE  METHOD 


4-  D  (SCRIPT  I  VS  NOTSS  (Tyyt  •!  n^m*  m4  AmAmIm 

Manuscript  for  Publication  (Technical  Report) 


a.  AUTMOAf4>  (L—t  MM,  Ami  mm,  mm 

Bartels,  Richard  H.  and  Golub,  Gene  H. 


4.  At  PORT  OATC 

June  2,  1967 


l«.  CONTRACT  ON  4NANT  NO. 

Nonr-255(37) 

a  PNOJRCT  NO. 

Nr-044-211 


10.  A  V  A  IL  ARILITY/LIMITATION  MOTICkS 


TA  MR.  RP  RRPS 

12 


S«.  OAI4tNATOA'4  A  CROAT  NUMStRffl) 


Releasable  without  limitations  on  dissemination. 


IS,  SPONSORINO  IMUTARV  ACTIVITY 


Office  of  Naval  Research  Code  432 
Washington,  D.C.  20360 


I S.  ABSTRACT 


^  An  implementation, using  Gaussian  LU  decomposition  with  row  interchanges, 
of  Stiefel's  exchange  algorithm  for  determining  a  Chebyshev  solution  to  an 
overdetermined  system  of  linear  equations  is  presented.  The  implementation 
is  computationally  more  stable  than  those  usually  given  in  the  literature. 

A  generalization  of  Stiefel's  algorithm  is  developed  which  permits  the 
occasional  exchange  of  two  equations  simultaneously.  Finally,  some 
experimental  comparisons  are  offered.  I  ^  ^ - 


DD  .KK.  1473 


Unclassified _ 

Security  Classification 


Unclaitiil  lied 


link  c 


KIV  WORD* 


LINK  A 


NOLI  I  WT 


LINK  B 


NOLI  I  WT  I  NOLI  I  WT 


Exchange  Algorithm 
Overdetermined  Linear  Systems 
Chebyshev  Solutions 


INSTRUCTIONS 


1.  ORIOINATINO  ACTIVITY,  Enim  tha  1 

S^SSffSBttSSSBBuSi  suthor)  L.ulng 

*to  BECUHTY  CLASSIFICATION:  Entar  tha  over- 

^.VlrtSS ©tottliwiSSl  K?; i»  «coU 

•new  with  epproprl.te  security  regulations. 

Ent*^ 

saa?Kfs..“Si  *  w  •  -  «■** 4  ■■  •",h~ 

“'report 

c.pu.1  tolly  1 }«•■ '  ‘ ,  UV.ucud  Without  Classifies- 

SA.’SSTffll  “ .JSE*.  •" M 

Immediately  following  the  UtU. 

a  DESCRIPTIVE  NOTES:  If  appropriate,  enter  the  type  ol 

'irvV‘^'f;ciu“rvfd.v.?wVe'n  yssSuZZm*  *• 

AUTHOR(S):  Enter  the  neme(e)  of  euthorf*) 

■  J  ShiT  Enter  last  name,  firet  name,  middle  initial. 

.  r  in  the  report.  Enter  taei  name.  Th*  name  of 

——  »»*--• 

(  on  the  report,  uee  date  of  publication. 

!  ,  ~Tti  NUMBER  OF  PAGES:  The  total  page  count 

R----  — ,h* 

number  of  pagea  containing  Information. 

NUMBER  OF  REFERENCES:  Enter  the  total  number  of 
I o-ferencee  cited  In  the  report. 

u  r  nNTRACT  OR  GRANT  NUMBER:  If  appropriate,  enter 

f£  .i  .<•  «»"““> «  *h‘ch 

the  report  waa  written. 

„  ...  oprtiECT  NUMBER:  Enter  the  appropriate 

...  ORIGINATOR’S 

a  ass  —■ . . 

be  unique  to  thla  report. 

-i 

«>■ •-««  ,hu  number(u)' 

l0.  AVAILABILITY/LIMITATION  NOTICES:  Enter  -» ^ 

ltatlonu  on  further  dissemination  of  the  report,  o.nc 


a  a  awnw 

Impoaed  by  aecurlty  claaalflcation,  uaing  standard  atatementa 
auch  aa: 

(1)  "Qualified  requeetere  may  obtain  coplaa  of  thla 
report  from  DDC.” 

(2)  "Foreign  announcement  and  dlaaamlnatlon  of  thla 
report  by  DDC  la  not  authorised.  ” 

(3)  “U.  S.  Government  agencloa  may  obtain  cop  lee  of 
thla  report  directly  from  DDC.  Other  qualified  DDC 
uaara  ehall  requeat  through 

•• 

_ a 

(41  "U.  S.  military  aganclea  may  obtain  coplaa  of  thla 

report  directly  from  DDC  Other  qualified  ueera 
ahell  requeat  through 


(5)  "Ail  dlatrlbutlon  of  thla  raport  la  controlled.  f>»al* 

Iflad  DDC  uaara  ahall  requaat  through 

m 

_ a 

If  tha  report  hee  been  furnlahed  to  the  Office  of  Technical 
Services,  Department  of  Commerce,  for  eale  to  the  public,  Indi¬ 
cate  thla  fact  end  enter  the  price,  If  known. 

IL  SUPPLEMENTARY  NOTES:  Uae  for  additional  explana¬ 
tory  notes. 

12.  SPONSORING  MILITARY  ACTIVITY:  Erter  the  name  of 
the  departmental  project  office  or  laboratory  aponaorlng  (par 
inf  for)  (h«  rt»e«rch  *nd  development  Include  eddreet* 

13  ABSTRACT:  Enter  an  abstract  giving  a  brief  and  factual 
summary  of  the  document  indicative  of  the  report,  even  though 
it  may  oleo  appear  elaewhere  in  the  body  of  the  technical  re- 
port.  If  additional  spaca  la  required,  a  continuation  ahaet  ahall 
be  attached. 

It  ie  highly  dealreble  that  the  abstract  of  classified  reports 
be  unclassified.  Each  paragraph  of  the  abstract  ehall  and  with 
an  indication  of  tha  military  sacurity  classification  of  tha  in¬ 
formation  in  tha  paragraph,  represented  aa  (T3).  <S),  (C),  er  (U). 

There  ie  no  limitation  on  the  length  of  the  abetract.  How¬ 
ever,  the  euggeated  length  is  from  150  to  225  words, 

14.  KEY  WORDS:  Key  words  are  technically  meaningful  tanas 
or  short  phraaea  that  characterize  e  report  end  may  be  used  aa 
index  entries  for  cataloging  the  report.  Key  words  must  bo 
selected  eo  that  no  aecurlty  claaalflcation  la  required.  Identi¬ 
fiers,  auch  aa  equipment  model  designation,  trade  name,  military 
project  code  name,  geographic  location,  may  be  used  as  key 
words  but  will  be  followed  by  an  indication  of  technical  con¬ 
text.  The  assignment  of  links,  roles,  and  weights  la  optional. 


RECEIVED 

SEP  2  9 1967 
CFSTI 


Erratum:  Insert  in  §8  before  the  Algol  6o  procedure 
of  Computer  Science  Report  No.  6?,  Stanford 
University. 

The  parameters  to  procedure  Chebyshev  are: 


identifier 

type 

comments 

m 

integer 

Number  of  equations . 

n 

integer 

Number  of  unknowns. 

A 

real  array 

Matrix  of  coefficients. 

Array  bounds  -  [0:m-l,  0:n-l]. 

d 

real  array 

Right-hand-side  vector. 

Array  bounds  -  [0:m-l]. 

h 

real  array 

Solution  vector. 

Array  bounds  -  [0:n-l]. 

ref set 

integer  array 

Final  reference  equation  numbers. 
Array  bounds  -  [0:n]. 

epz 

real 

Final  reference  deviation. 

zerolambda 

label 

Exit  for  condition  1  failure. 

insuf f ic ientrank 

label 

Exit  for  condition  2  failure, 
or  in  case  rank(A)  <  n  . 

The  parameters  m,  n,  A,  and  d  are  not  changed  by  Chebyshev. 
We  direct  the  user's  attention  to  the  identifier  eta  appearing  in  the 
procedure  and  to  the  comment  exp] -lining  its  value  and  purpose. 


